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Abstract 

The  desire  to  operate  denser  networks  in  order  to  monitor  seismic  activity  at  lower  thresholds 
leads  to  greater  emphasis  on  automated  data  processing.  An  algorithm  for  detecting  and 
characterizing  long-period  Rayleigh  wave  arrivals  has  been  developed  and  tested.  The  routine 
continuously  monitors  all  directions  of  approach  to  a  station,  in  a  manner  similar  to  beamforming. 
The  detector  is  based  on  cross-powers  between  the  Hilbert-transformed  vertical  and  rotated 
horizontal  signals,  so  it  is  sensitive  to  both  the  power  and  polarization  properties  of  the  three- 
component  wavefield.  Elliptically  polarized  Rayleigh  arrivals  are  enhanced,  while  linearly 
polarized  Love  waves  and  body  phases  are  suppressed.  A  test  using  one  month  of  data  from 
station  ANMO  demonstrated  that  this  technique  can  with  high  reliability  detect  Rayleigh  arrivals 
which  would  be  visible  to  an  analyst.  The  measured  arrival  times  and  azimuths  are  accurate 
enough  to  permit  automated  association  of  the  detections  to  events  in  a  bulletin. 
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Objective: 


Effective  verification  of  a  Comprehensive  Test  Ban  will  require  an  extensive  seismic  network 
capable  of  monitoring  at  a  fairly  low  magnitude  threshold.  Such  a  network  will  generate  a  large 
quantity  of  waveform  data  and  record  hundreds  or  even  thousands  of  events  daily.  The  resulting 
data  analysis  burden  will  force  an  increased  reliance  on  automated  processing.  Current  networks 
routinely  process  short-period  (SP)  waveforms  by  computer  to  automatically  detect  and  associate 
body  arrivals,  and  then  locate  the  events.  Improvements  to  these  routines  are  necessary  in  order  to 
improve  their  reliability.  In  contrast,  the  analysis  of  long-period  (LP)  waveforms,  important  for 
event  identification  via  the  M^mj,  discriminant,  remains  a  largely  manual  process.  First-generation 
algorithms  for  automated  LP  processing  have  only  recently  become  available.  One  approach  to 
this  problem,  currently  being  tested  by  the  prototype  International  Data  Center  under  GSETT-3, 
is  to  examine  the  LP  signals  for  the  existence  of  surface  waves  from  known  events.  An  alternate 
approach,  which  has  been  developed  and  tested  for  this  report,  is  to  continuously  scan  the  LP  data 
for  surface  arrivals  from  anywhere.  Detected  arrivals  can  then  be  associated  to  existing  event 
solutions,  or  used  to  help  constrain  the  solutions  for  events  which  produce  few  SP  detections. 

The  Rayleigh  detection  algorithm  discussed  here  exploits  the  elliptical  polarization  of  Rayleigh 
arrivals  both  to  enhance  detection  sensitivity  over  that  achieved  from  vertical  signals  alone,  and  to 
decrease  the  sensitivity  to  unwanted  body  or  Love  arrivals.  Several  authors  have  developed 
algorithms  for  polarization  analysis  of  seismic  waves.  Some  of  these  ( e.g .,  Simons,  19S8;  Samson, 
1977)  have  required  knowledge  of  event  location,  so  that  the  radial  and  transverse  directions  are 
known  a  priori.  Others  {e.g.,  Samson  and  Olson,  1981)  have  relied  on  continuously  tracking  the 
dominant  polarization  direction,  in  essence  dynamically  rotating  the  waveforms  to  this  continually 
changing  azimuth.  The  method  addressed  in  this  work  seeks  to  continuously  and  independently 
monitor  all  directions  for  Rayleigh  arrivals.  The  technique  is  comparable  to  beamforming  done  for 
seismic  arrays,  in  which  beams  are  pointed  in  several  different  directions  then  separately  scanned 
for  arrivals.  Simultaneously  watching  all  directions  should  offer  signal-to-noise  improvement 
compared  to  tracking  only  the  current  dominant  direction. 


Research  Accomplished; 

An  automated  Rayleigh  wave  detection  routine  has  been  designed,  and  its  performance  tested 
using  one  month  of  LP  data.  Three-component  long-period  or  broadband  waveforms  are  required 
for  the  Rayleigh  detection  algorithm.  The  first  step  in  the  processing  is  to  detrend  and  then 
bandpass  filter  the  signals.  A  third-order  Butterworth  filter  with  a  passband  of  0.01  to  0.06  Hz  has 
produced  good  results.  Next,  the  vertical  signal  is  Hilbert  transformed  using  a  255-point  FIR  filter 
(Steams  and  David,  1988).  This  has  the  effect  of  phase  shifting  the  vertical  waveform  by  90°, 
which  converts  the  elliptical  polarization  of  Rayleigh  waves  into  linear  motion.  A  beneficial  side 
effect  of  this  operation  is  to  change  linearly  polarized  body  waves  to  elliptical  polarization,  which 
will  de-emphasize  them  in  subsequent  processing.  At  each  time  step,  the  3-by-3  data  covariance 
matrix  is  calculated  over  a  specified  window  length.  A  duration  of  80  s  has  been  used  for  this 
short-term-average  (ST A)  window,  equivalent  to  4  cycles  at  the  dominant  Rayleigh  period  near 
20  s.  The  covariance  matrix  at  each  time  step  is  then  rotated  over  azimuth  from  0°  to  85°  in  5° 
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increments.  At  each  azimuth  step,  the  two  covariance  elements  which  represent  the  cross-powers 
between  the  vertical  and  each  of  the  rotated  horizontals  are  extracted.  These  are  used  as  the  STA 
estimates  for  the  two  orthogonal  azimuths  corresponding  to  the  angle  of  rotation  (0°  to  85°)  and 
the  direction  90°  clockwise  from  this  (90°  to  175°).  It  is  not  necessary  to  calculate  the 
covariances  for  azimuths  of  180°  to  355°,  since  the  values  would  simply  be  the  negative  of  those 
for  0°  to  175°.  A  strong  positive  STA  response  indicates  Rayleigh  energy  with  the  appropriate 
retrograde  elliptical  motion  arriving  from  the  specified  direction.  Using  the  cross-power  between 
the  vertical  and  horizontal  motions  should  improve  the  signal-to-noise  by  about  a  factor  of  two 
(3  dB),  compared  to  using  the  power  on  either  axis  alone.  Long-term-averages  (LTAs)  of  each  of 
the  ST  As  are  updated  at  each  time  step.  A  recursive  scheme  was  used  for  the  LTA  computations, 
providing  exponential  weighting  of  past  ST  As.  The  time  constant  for  this  weighting  was  set  at 
800  s. 

A  detection  function  is  produced  for  each  look  direction  by  calculating  a  weighted  STA/LTA 
ratio  versus  time.  Alternatively,  one  can  view  the  result  as  a  detection  surface  with  an  altitude  that 
varies  as  a  function  of  both  azimuth  and  time.  Two  weights  are  applied  to  the  STA/LTA  ratios  of 
the  cross-powers.  The  first  of  these  weights  is  the  absolute  value  of  the  correlation  coefficient 
between  the  vertical  and  the  radial  signal  for  a  given  azimuth.  Because  the  vertical  has  been 
Hilbert  transformed,  it  should  be  well  correlated  with  the  horizontal  when  a  Rayleigh  wave  is 
arriving  from  the  specified  direction.  The  second  weight  is  based  on  the  relative  amplitudes  of  the 
vertical  and  radial  components.  Its  purpose  is  to  emphasize  arrivals  whose  Z  and  R  motions  are 
similar  in  amplitude.  Both  of  these  weights  are  readily  calculated  from  the  elements  of  the  rotated 
covariance  matrix.  The  detection  surface  is  scanned  for  peaks  exceeding  a  threshold  altitude. 
When  a  detection  is  made,  the  peak's  onset  time,  amplitude  and  central  azimuth  are  recorded, 
along  with  the  peak  amplitude  of  the  vertical  signal  following  the  detection.  The  entire  Rayleigh 
detection  algorithm  was  implemented  in  MATLAB,  so  it  can  be  easily  modified  to  test  different 
processing  details.  The  current  routine  should  be  able  to  process  full-period  LP  data  (1  sps)  from 
roughly  100  stations  on  a  single  Pentium-class  PC. 

Figure  1(a)  shows  20  min  of  LP  data  from  station  COL  near  Fairbanks,  AK  for  an  earthquake  in 
the  Aleutians.  The  horizontal  components  have  been  rotated  to  the  radial  (middle  trace)  and 
transverse  (bottom  trace)  directions.  A  strong  Rayleigh  wave  arrives  near  600  s  on  the  LPZ  and 
LPR  traces.  The  LPT  signal  shows  a  strong  Love  arrival  near  400  s,  with  coda  extending  into  the 
Rayleigh  window.  This  coda  contaminates  a  standard  polarization  analysis  of  the  Rayleigh  wave. 
Figure  1(b)  displays  the  detection  surface  produced  by  these  records;  the  surface  is  shown  for 
azimuths  of  1 80°  to  360°,  where  the  data  produce  a  positive  peak.  The  Rayleigh  wave  produced  a 
very  strong  peak  at  the  proper  time  and  azimuth  for  this  event.  Note  that  the  large  Love  arrival 
had  little  effect  on  the  surface.  The  Love  coda  caused  the  ridge  running  to  the  left  from  the  peak, 
but  did  not  significantly  bias  the  peak's  location.  Figure  2(a)  shows  a  much  weaker  Rayleigh 
arrival  at  station  ANMO  (Albuquerque,  NM)  from  an  event  in  the  South  Pacific.  For  this  plot  the 
traces  were  bandpass  filtered  between  0.01  Hz  and  0.06  Hz  to  enhance  the  arrival.  Again  the 
horizontals  have  been  rotated  to  radial  and  transverse,  and  the  Rayleigh  wave  begins  just  prior  to 
600  s.  The  resulting  detection  surface  in  Figure  2(b)  shows  a  prominent  peak  at  the  Rayleigh 
onset  time,  centered  near  the  correct  backazimuth  to  the  event.  This  peak  jumps  abruptly  from  the 
preceding  terrain,  despite  the  relatively  low  SNR  of  the  arrival  on  the  seismograms.  A  second 
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peak  occurs  near  900  s  on  the  plot,  caused  by  the  resurgence  in  the  Rayleigh  coda  which  can  be 
seen  on  the  LPZ  trace. 

To  test  the  performance  of  the  detection  algorithm  in  routine  operation,  I  processed  all  of  the 
ANMO  data  for  the  month  of  October,  1993.  The  resulting  detection  log  was  compared  to  the 
NEIC  event  bulletin  for  this  month.  A  total  of  192  detections  were  reported,  of  which  157  could 
be  attributed  to  108  different  events  in  the  bulletin.  For  those  events  which  caused  multiple 
detections,  the  trigger  corresponding  to  the  Rayleigh  onset  could  usually  be  readily  determined 
from  its  time  and  azimuth.  Multiple  detections  for  an  event  were  usually  due  to  fluctuations  in  the 
Rayleigh  coda,  overtone  energy  preceding  the  Rayleigh,  or  S  arrivals.  Errors  in  the  estimated 
backazimuths  for  identified  Rayleigh  arrivals  had  a  mean  of  0°  and  a  standard  deviation  of  15°. 
Thus  the  backazimuth  measurements  are  sufficiently  accurate  to  assist  automated  association  of 
Rayleigh  detections  with  events.  Figure  3  shows  the  performance  of  the  detector  as  a  function  of 
event  magnitude  and  distance,  for  all  events  which  had  Ms  values  in  the  NEIC  bulletin.  Solid 
circles  represent  events  which  produced  Rayleigh  detections  by  the  algorithm,  open  circles  events 
which  did  not.  This  plot  indicates  that  the  routine  can  detect  Rayleigh  waves  from  most  events 
with  Ms  over  4  to  distances  of  30°,  with  the  sensitivity  gradually  decreasing  to  an  Ms  of  about  5 
at  the  flirthest  ranges.  The  undetected  large  events  at  a  distance  of  108°  were  aftershocks  in  the 
coda  of  larger  events.  Adjusting  the  LTA  time  constant  or  changing  the  logic  for  ending  a 
previous  detection  may  yield  better  performance  in  such  cases. 


Conclusions,  and  Recommendations: 

An  algorithm  for  detecting  Rayleigh  surface  wave  arrivals  on  three-component  long-period 
records  has  been  developed  and  tested.  The  routine  continuously  monitors  all  directions  of 
approach  to  a  station,  in  a  manner  analogous  to  array  beamforming.  The  detector's  response  is 
sensitive  to  both  the  power  and  polarization  properties  of  the  wavefield,  so  that  Rayleigh  waves 
are  preferentially  enhanced  over  other  arrivals.  A  one-month  test  of  the  algorithm  with  station 
ANMO  demonstrated  that  it  can  reliably  detect  an  overwhelming  majority  of  the  Rayleigh  waves 
visible  on  the  seismograms,  and  the  detection  time  and  azimuth  measures  are  sufficiently  accurate 
for  reliable  association  to  bulletin  events. 

Further  work  is  needed  in  three  areas.  First,  estimates  of  arrival  amplitude  and  period  need  to  be 
refined,  in  order  to  enable  accurate  automated  Mg  measurements.  Second,  better  methods  for 
screening  non-Rayleigh  detections  should  be  developed  and  tested.  Examining  a  detection's 
frequency  content  may  be  useful  here.  Finally,  more  testing  is  needed  to  establish  the  performance 
for  regional  events  (A  <  30°),  which  will  become  of  increasing  importance  under  a  CTBT. 


References: 

Samson,  J.  C.  (1977).  Matrix  and  Stokes  vector  representations  of  detectors  for  polarized 
waveforms:  theory,  with  some  applications  to  teleseismic  waves,  Geophys.  J.  R.  astr.  Soc. 
51,583-603. 


853 


Samson,  J.  C.  and  J.  V.  Olson  (1981).  Data-adaptive  polarization  filters  for  multichannel 
geophysical  data.  Geophysics  46,  1423-1431. 

Simons,  R.  S.  (1968).  A  surface  wave  particle  motion  discrimination  process.  Bull.  Seism.  Soc. 
Am.  58,  629-637. 

Steams,  S.  D.  and  R.  A.  David  (1988).  Signal  Processing  Algorithms ,  Prentice-Hall,  Inc., 
Englewood  Cliffs,  NJ. 


854 


Detector  Output  for  COL,  5/8/86 
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Figure  1 .  (a)  Vertical,  radial  and  transverse  records  from  station  COL  of  an  Mg  5.4  earthquake  in  the  Aleutian  Islands  on 
May  8,  1986,  at  a  distance  of  20°  and  a  backazimuth  of  241°.  (b)  Output  of  the  detection  algorithm  for  the  signals 
shown  in  (a). 
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Figure  2.  (a)  Vertical,  radial  and  transverse  records  from  station  ANMO  of  an  Ms  4.3  earthquake  in  the  South  Pacific  on 
October  7,  1993,  at  a  distance  of  97°  and  a  backazimuth  of  232°.  (b)  Output  of  the  detection  algorithm  for  the  signals 
shown  in  (a). 
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Figure  3.  Detector  performance  as  a  function  of  event  distance  and  surface  wave  magnitude.  Solid 
circles  represent  events  whose  Rayleigh  waves  were  detected  by  the  algorithm,  open  circles 
indicate  that  a  Rayleigh  arrival  was  not  detected. 
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